C-------------------------------------------------------------------------------
C  x: Mx=y
C-------------------------------------------------------------------------------
      subroutine precond(prec, Nx, Ny, Nz, A, D, y, x)
      implicit none
      
      integer Nx, Ny, Nz, prec
      
      double precision A(Nx*Ny*Nz,4), D(Nx*Ny*Nz), y(Nx*Ny*Nz)
      double precision x(Nx*Ny*Nz)
      
      double precision, dimension(:), allocatable :: q
      
      integer NxNy, Nt, i
      
      allocate(q(Nx*Ny*Nz))
      
      NxNy = Nx*Ny
      Nt   = Nx*Ny*Nz
      

      if (prec .eq. 1) then
C--------------------
C  Diagonal scaling
C--------------------
      	do i=1,Nt
        	x(i) = y(i) / A(i,1)
      	end do
      	return
      elseif (prec .eq. 2) then
C-----------------------
C  Incomplete Choleski
C-----------------------
      do i=1,1
	q(i) = y(i) / D(i)
      end do
      do i=2,Nx
	q(i) = (y(i) - A(i-1,2)*q(i-1)) / D(i)
      end do
      do i=Nx+1,NxNy
	q(i) = (y(i) - A(i-1,2)*q(i-1) - A(i-Nx,3)*q(i-Nx)) / D(i)
      end do
      do i=NxNy+1,Nt
	q(i) = (y(i) - A(i-1,2)*q(i-1) - A(i-Nx,3)*q(i-Nx)
     1                                 - A(i-NxNy,4)*q(i-NxNy)) / D(i)
      end do
      
      do i=Nt,Nt,-1
        x(i) = q(i)
      end do
      do i=Nt-1,Nt-Nx+1,-1
        x(i) = q(i) - A(i,2)*x(i+1) / D(i)
      end do
      do i=Nt-Nx,Nt-NxNy+1,-1
        x(i) = q(i) - (A(i,2)*x(i+1) + A(i,3)*x(i+Nx)) / D(i)
      end do
      do i=Nt-NxNy,1,-1
        x(i) = q(i) - (A(i,2)*x(i+1) + A(i,3)*x(i+Nx)
     1                               + A(i,4)*x(i+NxNy)) / D(i)
      end do
      
C      write(6,"(E25.15)") (x(i), i=1,Nt)
  
      deallocate(q)
      
      return
      endif

      end
C-------------------------------------------------------------------------------
C  
C-------------------------------------------------------------------------------
      subroutine make_cholinc(Nx, Ny, Nz, A, D)
      implicit none
      
      integer Nx, Ny, Nz
      
      double precision A(Nx*Ny*Nz,4)
      double precision D(Nx*Ny*Nz)
      
      integer NxNy, Nt, i
      
      NxNy = Nx*Ny
      Nt   = Nx*Ny*Nz
      
      do i=1,1
        D(i) = A(i,1)
      end do
      
      do i=2,Nx
        D(i) = A(i,1) - A(i-1,2)**2/D(i-1)
      end do
      
      do i=Nx+1,NxNy
        D(i) = A(i,1) - A(i-1,2)**2/D(i-1) - A(i-Nx,3)**2/D(i-Nx)
      end do
      
      do i=NxNy+1,Nt
        D(i) = A(i,1) - A(i-1,2)**2/D(i-1) - A(i-Nx,3)**2/D(i-Nx)
     1                                   - A(i-NxNy,4)**2/D(i-NxNy)
      end do
      
C      write(6,"(E25.15)") (D(i), i=1,Nx*Ny*Nz)
  
      return
      end
